Diabetes
2023-03-01
library(readxl)
mydata <- read_xlsx("./Diabetes.xlsx")
mydata <- as.data.frame(mydata)
head(mydata)## ID Pregnancies Glucose Pressure BMI Age Diabetes
## 1 1 1 97 64 18.2 21 neg
## 2 2 1 83 68 18.2 27 neg
## 3 3 1 97 70 18.2 21 neg
## 4 4 0 104 76 18.4 27 neg
## 5 5 1 80 55 19.1 21 neg
## 6 6 3 99 80 19.3 30 negDescription of variables used:
- ID: Patient ID.
- Pregnancies: How many times the women was Pregnancies.
- Glucose: Level of glucose in mg/dl.
- Pressure: Diastolic blood pressure in mm HG.
- BMI: Body BMI index (kg/(m2)).
- Age in years.
- Diabetes: Test for diabetes (positive or negative)
mydata$Diabetes <- factor(mydata$Diabetes,
levels = c("neg", "pos"),
labels = c("neg", "pos"))library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionmydata %>%
select(Diabetes, Pregnancies, Glucose, Pressure, BMI, Age) %>%
summary()## Diabetes Pregnancies Glucose Pressure
## neg:422 Min. :0.000 Min. : 44.0 Min. : 24.0
## pos:184 1st Qu.:1.000 1st Qu.: 99.0 1st Qu.: 64.0
## Median :2.000 Median :115.0 Median : 70.0
## Mean :2.713 Mean :120.7 Mean : 71.4
## 3rd Qu.:4.000 3rd Qu.:139.0 3rd Qu.: 80.0
## Max. :7.000 Max. :199.0 Max. :122.0
## BMI Age
## Min. :18.20 Min. :21.00
## 1st Qu.:27.40 1st Qu.:23.00
## Median :32.00 Median :27.00
## Mean :32.39 Mean :31.12
## 3rd Qu.:36.58 3rd Qu.:36.00
## Max. :67.10 Max. :70.00fit0 <- glm(Diabetes ~ 1, #Dependent and explanatory variables
family = binomial, #Binary logistic regression
data = mydata)
summary(fit0)##
## Call:
## glm(formula = Diabetes ~ 1, family = binomial, data = mydata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.83007 0.08834 -9.396 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 744.06 on 605 degrees of freedom
## Residual deviance: 744.06 on 605 degrees of freedom
## AIC: 746.06
##
## Number of Fisher Scoring iterations: 4exp(cbind(odds = fit0$coefficients, confint.default(fit0))) #Odds for Y=1## odds 2.5 % 97.5 %
## (Intercept) 0.436019 0.3666974 0.5184452head(fitted(fit0)) #Estimated probability for Y=1## 1 2 3 4 5 6
## 0.3036304 0.3036304 0.3036304 0.3036304 0.3036304 0.3036304ACC <- 422/606 #Prop. of correctly classified units
ACC## [1] 0.6963696fit1 <- glm(Diabetes ~ Pregnancies + Glucose + Pressure + BMI + Age,
family = binomial,
data = mydata)library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodevif(fit1)## Pregnancies Glucose Pressure BMI Age
## 1.251598 1.044367 1.202577 1.136284 1.370504mydata$StdResid <- rstandard(fit1)
hist(mydata$StdResid,
main = "Histogram of standardised residuals",
ylab = "Frequency",
xlab = "Standardised residuals")
library(dplyr)
# Calculate Cook's distance
mydata$CooksD <- cooks.distance(fit1)
# Define Cook's distance threshold
cd_threshold <- 4/nrow(mydata)
# Display apartments with Cook's distance above the threshold
influential_units <- mydata %>%
filter(CooksD > cd_threshold) %>%
arrange(desc(CooksD)) %>%
select(ID, Diabetes, Pregnancies, Glucose, Pressure, BMI, Age, CooksD)
head(influential_units, 6) %>%
mutate(CooksD = round(CooksD, 3))## ID Diabetes Pregnancies Glucose Pressure BMI Age CooksD
## 1 590 neg 0 173 78 46.5 58 0.031
## 2 173 pos 0 105 84 27.9 62 0.023
## 3 30 neg 0 161 50 21.9 65 0.021
## 4 600 neg 0 165 90 52.3 23 0.016
## 5 382 neg 7 179 95 34.2 60 0.016
## 6 603 pos 1 88 30 55.0 26 0.016fit1 <- glm(Diabetes ~ Pregnancies + Glucose + Pressure + BMI + Age,
family = binomial,
data = mydata)anova(fit0, fit1, test = "Chisq") #Comparision of models based on -2LL statistics## Analysis of Deviance Table
##
## Model 1: Diabetes ~ 1
## Model 2: Diabetes ~ Pregnancies + Glucose + Pressure + BMI + Age
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 605 744.06
## 2 600 551.24 5 192.82 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1logLik_fit0 <- as.numeric(logLik(fit0))
logLik_fit1 <- as.numeric(logLik(fit1))
McFadden_R2 <- 1 - (logLik_fit1 / logLik_fit0)
McFadden_R2## [1] 0.2591475summary(fit1)##
## Call:
## glm(formula = Diabetes ~ Pregnancies + Glucose + Pressure + BMI +
## Age, family = binomial, data = mydata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.580862 0.870336 -9.859 < 2e-16 ***
## Pregnancies 0.066884 0.053407 1.252 0.21045
## Glucose 0.034899 0.003975 8.779 < 2e-16 ***
## Pressure -0.013465 0.009253 -1.455 0.14560
## BMI 0.094391 0.017034 5.541 3e-08 ***
## Age 0.031009 0.010846 2.859 0.00425 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 744.06 on 605 degrees of freedom
## Residual deviance: 551.24 on 600 degrees of freedom
## AIC: 563.24
##
## Number of Fisher Scoring iterations: 5exp(cbind(OR = fit1$coefficients, confint.default(fit1))) #Odds ratio for Y=1 (with 95% CI)## OR 2.5 % 97.5 %
## (Intercept) 0.0001876632 3.408375e-05 0.001033263
## Pregnancies 1.0691712041 9.629124e-01 1.187155800
## Glucose 1.0355147782 1.027478e+00 1.043614584
## Pressure 0.9866249190 9.688935e-01 1.004680871
## BMI 1.0989891661 1.062905e+00 1.136298837
## Age 1.0314952587 1.009800e+00 1.053656253mydata$EstProb <- fitted(fit1) #Estimates of probabilities for Y=1: P(Y=1)
head(mydata[c("ID", "Diabetes", "EstProb")])## ID Diabetes EstProb
## 1 1 neg 0.02604660
## 2 2 neg 0.01838142
## 3 3 neg 0.02407367
## 4 4 neg 0.03227302
## 5 5 neg 0.01783469
## 6 6 neg 0.03731102mydata$Classification <- ifelse(test = mydata$EstProb < 0.50,
yes = "NO",
no = "YES")
#If estimated probability is below 0.50, person is classified into a group of people without problems of loan repayment (NO), otherwise YES.
mydata$Classification <- factor(mydata$Classification,
levels = c("NO", "YES"),
labels = c("NO", "YES"))Table <- table(
Actual = mydata$Diabetes,
Predicted = mydata$Classification
)
Table## Predicted
## Actual NO YES
## neg 382 40
## pos 89 95# Install once if not already installed
# install.packages("pROC")
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var# Loads the pROC package, which provides functions to compute and plot
# Receiver Operating Characteristic (ROC) curves and to calculate AUC.
# Compute ROC curve
roc_fit1 <- roc(response = mydata$Diabetes,
predictor = mydata$EstProb,
levels = c("neg", "pos"),
direction = "<")
# Constructs the ROC curve object for model fit1.
# - response: the observed outcome (diabetes status),
# - predictor: predicted probabilities P(Y = 1) from the logistic model,
# - levels: specifies which category is treated as non-event ("neg")
# and which as event ("pos"),
# - direction: indicates that higher predicted probabilities correspond
# to a higher likelihood of the event (diabetes).# Area Under the Curve (AUC)
auc_fit1 <- auc(roc_fit1)
auc_fit1## Area under the curve: 0.8336# Computes the area under the ROC curve.
# The AUC summarises the model’s overall discriminatory ability
# independently of any particular classification threshold.
# Plotting ROC curve, including AUC
plot(roc_fit1,
col = "blue",
lwd = 2,
main = paste("ROC Curve for Model fit1 (AUC =",
round(auc_fit1, 3), ")"))
# Plots the ROC curve.
# - The y-axis represents sensitivity (true positive rate),
# - The x-axis represents 1 - specificity (false positive rate),
# - The AUC value is displayed in the plot title for interpretation.
# Add reference line for random classification
abline(a = 0, b = 1, lty = 2, col = "gray")
# Adds the 45-degree reference line corresponding to random classification.
# A model with no discriminatory power would lie close to this line.coords(roc_fit1,
x = "best",
best.method = "youden",
ret = c("threshold", "sensitivity", "specificity", "accuracy"))## threshold sensitivity specificity accuracy
## threshold 0.2382797 0.8152174 0.7014218 0.7359736